{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "38d84418",
   "metadata": {},
   "source": [
    "# Simulating Single Phase Transport \n",
    "\n",
    "The point of an OpenPNM simulation is ultimately to compute some transport process.  In this notebook we will cover the following subjects:\n",
    "\n",
    "- Defining conductance \n",
    "- Settings boundary conditions \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac063b0",
   "metadata": {},
   "source": [
    "Start by defining a network.  We'll use the ``Demo`` class which happens to include all the geometrical pore-scale models already. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "94917bac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "net : <openpnm.network.Demo at 0x25645ed9b30>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.coordination_number                                          25 / 25\n",
      "  2  pore.coords                                                       25 / 25\n",
      "  3  pore.diameter                                                     25 / 25\n",
      "  4  pore.max_size                                                     25 / 25\n",
      "  5  pore.seed                                                         25 / 25\n",
      "  6  pore.volume                                                       25 / 25\n",
      "  7  throat.conns                                                      40 / 40\n",
      "  8  throat.cross_sectional_area                                       40 / 40\n",
      "  9  throat.diameter                                                   40 / 40\n",
      " 10  throat.diffusive_size_factors                                     40 / 40\n",
      " 11  throat.hydraulic_size_factors                                     40 / 40\n",
      " 12  throat.length                                                     40 / 40\n",
      " 13  throat.lens_volume                                                40 / 40\n",
      " 14  throat.max_size                                                   40 / 40\n",
      " 15  throat.spacing                                                    40 / 40\n",
      " 16  throat.total_volume                                               40 / 40\n",
      " 17  throat.volume                                                     40 / 40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.back                                                               5\n",
      "  2  pore.domain_1                                                          25\n",
      "  3  pore.front                                                              5\n",
      "  4  pore.left                                                               5\n",
      "  5  pore.right                                                              5\n",
      "  6  pore.surface                                                           16\n",
      "  7  throat.domain_1                                                        40\n",
      "  8  throat.surface                                                         16\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import openpnm as op\n",
    "op.visualization.set_mpl_style()\n",
    "\n",
    "pn = op.network.Demo(shape=[5, 5, 1], spacing=5e-5)\n",
    "print(pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7a0902f",
   "metadata": {},
   "source": [
    "## Define Phase Viscosity\n",
    "\n",
    "To fully illustrate the process of performing transport calculations, we'll use an empty ``Phase`` object and add all the needed properties manually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b5e3a75a",
   "metadata": {},
   "outputs": [],
   "source": [
    "water = op.phase.Phase(network=pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8973bf1",
   "metadata": {},
   "source": [
    "Let's assume that we are interested in pressure driven flow. This requires knowing the viscosity of the phase, so let's add a pore-scale model for computing the viscosity of water:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "029d1797",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "phase_01 : <openpnm.phase.Phase at 0x2564c657180>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.pressure                                                     25 / 25\n",
      "  2  pore.temperature                                                  25 / 25\n",
      "  3  pore.viscosity                                                    25 / 25\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.all                                                               25\n",
      "  2  throat.all                                                             40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "water.add_model(propname='pore.viscosity',\n",
    "                model=op.models.phase.viscosity.water_correlation)\n",
    "print(water)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21aab8de",
   "metadata": {},
   "source": [
    "And we can check the individual values to verify they make sense:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8833d1ee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319\n",
      " 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319\n",
      " 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319\n",
      " 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319\n",
      " 0.00089319]\n"
     ]
    }
   ],
   "source": [
    "print(water['pore.viscosity'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95ba0972",
   "metadata": {},
   "source": [
    "## Basic Conductance Calculation\n",
    "\n",
    "Determining the conductance of the conduits between each pair of connected pores is the most important part of performing a simulation. The details of conductance models are covered elsewhere. For this demonstration will assume the very simplest case where all pressure loss occurs in the throats.\n",
    "\n",
    "Recall the Hagan-Poiseuille equation for fluid flow through a cylindrical tube:\n",
    "\n",
    "$$ Q = \\frac{\\pi R^4}{8 \\mu L} \\Delta P$$\n",
    "\n",
    "where $R$ and $L$ are the radius and length of the throat, and $\\mu$ is the viscosity of the fluid.  Together this prefactor can be referred to as the hydraulic conductance, $g_h$, giving:\n",
    "\n",
    "$$ Q = g_h \\Delta P $$\n",
    "\n",
    "So the aim is the compute values of $g_h$ for each throat.  We start by doing this manually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "647bde31",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00220815e-15 1.05133682e-15 2.65744048e-15 2.69017014e-15\n",
      " 9.59190878e-15 1.36708286e-14 6.66830487e-14 7.08737841e-14\n",
      " 1.99903177e-15 2.00595117e-15 4.79376978e-15 4.31484439e-15\n",
      " 3.19680625e-15 2.81277856e-15 5.16906031e-15 4.79857420e-15\n",
      " 8.55176371e-15 1.52083335e-15 1.23450963e-15 2.47324217e-15\n",
      " 9.86091431e-15 9.79143235e-16 3.16100872e-14 2.89691924e-15\n",
      " 3.71177597e-14 1.01109214e-14 1.90884769e-15 2.45494330e-14\n",
      " 5.38067275e-15 7.13185895e-15 2.56325956e-14 1.72969382e-15\n",
      " 2.28431461e-14 4.24865495e-15 5.99267859e-15 9.64675766e-15\n",
      " 2.86021000e-15 1.49648281e-15 2.35465588e-15 9.68639559e-15]\n"
     ]
    }
   ],
   "source": [
    "R = pn['throat.diameter']/2\n",
    "L = pn['throat.length']\n",
    "mu = water['throat.viscosity']  # See ProTip below\n",
    "water['throat.hydraulic_conductance'] = np.pi*R**4/(8*mu*L)\n",
    "print(water['throat.hydraulic_conductance'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7720d25b",
   "metadata": {},
   "source": [
    "```{admonition} Phase can do automatic interpolation to get throat values\n",
    "  The ``Phase`` class has a special ability to interpolate pore values to throats, and vice-versa. In PNM simulations all the balances are solved for each pore, so the thermodynamic properties like temperature, pressure, etc. are all defined on pores. Consequently, the physical properties are also defined in pores, like viscosity. However, as shown above we often want viscosity values in the throats. OpenPNM provides a shortcut for this, such that if you request a throat property that does not exist, it will attempt to fetch the pores values and do a linear interpolation of values to produce an array of throat values.  There is also a function for this, ``water.interpolate_data('throat.viscosity')``, but the ``water['throat.viscosity']`` shortcut is very convenient. The automatic interpolation can be disabled in the `phase.settings`.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d4544d8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water.interpolate_data('throat.viscosity')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cf30019f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,\n",
       "       0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water['throat.viscosity']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e71e725",
   "metadata": {},
   "source": [
    "## Create Algorithm Object\n",
    "\n",
    "OpenPNM contains a variety of ``Algorithm`` classes in the ``openpnm.algorithms`` module.  Let's initialize a ``StokesFlow`` algorithm, since this simulates pressure driven flow through the network. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "29b8bf9d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "stokes_01 : <openpnm.algorithms.StokesFlow at 0x2564c673f40>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.bc.rate                                                       0 / 25\n",
      "  2  pore.bc.value                                                      0 / 25\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.all                                                               25\n",
      "  2  throat.all                                                             40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "sf = op.algorithms.StokesFlow(network=pn, phase=water)\n",
    "print(sf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72c4d6c8",
   "metadata": {},
   "source": [
    "## Assign Boundary Conditions\n",
    "\n",
    "As can be seen in the print-out above there are predefined ``'pore.bc'`` arrays, but they contain no valid values, meaning they are all ``nans``. Once we set some boundary conditions, this will change.  Let's apply pressure BCs on one side of the network, and rate BCs on the other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3caad7c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "stokes_01 : <openpnm.algorithms.StokesFlow at 0x2564c673f40>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.bc.rate                                                       5 / 25\n",
      "  2  pore.bc.value                                                      5 / 25\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.all                                                               25\n",
      "  2  throat.all                                                             40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "sf.set_value_BC(pores=pn.pores('left'), values=100_000)\n",
    "sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)\n",
    "print(sf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e81f2e4b",
   "metadata": {},
   "source": [
    "```{admonition} All boundary conditions are preceeded with 'pore.bc'\n",
    "  All boundary conditions are stored as 'pore.bc.<type>', which means that OpenPNM's dictionary lookup tricks can be used to see all types and values of bc's using: `sf['pore.bc']` which will return a `dict`. This can be used as shown below:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f14cea28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['rate', 'value'])\n"
     ]
    }
   ],
   "source": [
    "print(sf['pore.bc'].keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7114b180",
   "metadata": {},
   "source": [
    "Now we can see there are 5 valid values of each type.  The ``sf`` algorithm will look for ``'throat.hydraulic_conductance'`` on ``water`` by default, so we can just ``run``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "732f855a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                       \r"
     ]
    }
   ],
   "source": [
    "soln = sf.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613df06f",
   "metadata": {},
   "source": [
    "The ``run`` method solves the mass balance around each pore and computes the pressure within each pore that is required to sustain the flow defined by the boundary conditions.  The ``soln`` object that is returned is a dictionary with the key corresponding to the quantity that was solved for.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cf47340a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'pore.pressure': SteadyStateSolution([100000.        , 100000.        , 100000.        ,\n",
      "                     100000.        , 100000.        , 111994.91786881,\n",
      "                     109021.80858086, 106007.58742924, 105515.18852917,\n",
      "                     104499.41168208, 126513.73814012, 120297.06671979,\n",
      "                     113401.97906307, 115761.92671225, 117822.18393838,\n",
      "                     132725.59341895, 133551.80327216, 120247.96545073,\n",
      "                     129309.20272742, 135160.97831439, 148957.44526658,\n",
      "                     155574.17873369, 168444.33285769, 161719.71531674,\n",
      "                     148786.89313432])}\n"
     ]
    }
   ],
   "source": [
    "print(soln)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd75fd50",
   "metadata": {},
   "source": [
    "The reason for the dict format is to provide a consistent API for single components and multiphysics, where multiple different quantities might be solved for.  However, these ``'pore.pressure'`` values are also written to the dictionary of the algorithm object as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "520819ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "stokes_01 : <openpnm.algorithms.StokesFlow at 0x2564c673f40>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.bc.rate                                                       5 / 25\n",
      "  2  pore.bc.value                                                      5 / 25\n",
      "  3  pore.initial_guess                                                25 / 25\n",
      "  4  pore.pressure                                                     25 / 25\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.all                                                               25\n",
      "  2  throat.all                                                             40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(sf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc617629",
   "metadata": {},
   "source": [
    "Finally we can look at how much pressure was required in the \"right\" pores to meet the required flow rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "16a3d0d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([148957.44526658, 155574.17873369, 168444.33285769, 161719.71531674,\n",
       "       148786.89313432])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sf['pore.pressure'][pn.pores('right')]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55cc13bb",
   "metadata": {},
   "source": [
    "So we can see that 150 kPa was required to accomplish the requested flow. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b9cba6",
   "metadata": {},
   "source": [
    "## Rigorous Conductance Calculation\n",
    "\n",
    "The above demonstration used a very simplistic conductance calculation.  It was also stated that computing conductance is the most important part of doing a PNM simulation.  To finish this notebook, we'll look more closely at this process."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "751a30eb",
   "metadata": {},
   "source": [
    "### Manual Method\n",
    "Let's print the network object again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "322b2746",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "══════════════════════════════════════════════════════════════════════════════\n",
      "net : <openpnm.network.Demo at 0x25645ed9b30>\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Properties                                                   Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.coordination_number                                          25 / 25\n",
      "  2  pore.coords                                                       25 / 25\n",
      "  3  pore.diameter                                                     25 / 25\n",
      "  4  pore.max_size                                                     25 / 25\n",
      "  5  pore.seed                                                         25 / 25\n",
      "  6  pore.volume                                                       25 / 25\n",
      "  7  throat.conns                                                      40 / 40\n",
      "  8  throat.cross_sectional_area                                       40 / 40\n",
      "  9  throat.diameter                                                   40 / 40\n",
      " 10  throat.diffusive_size_factors                                     40 / 40\n",
      " 11  throat.hydraulic_size_factors                                     40 / 40\n",
      " 12  throat.length                                                     40 / 40\n",
      " 13  throat.lens_volume                                                40 / 40\n",
      " 14  throat.max_size                                                   40 / 40\n",
      " 15  throat.spacing                                                    40 / 40\n",
      " 16  throat.total_volume                                               40 / 40\n",
      " 17  throat.volume                                                     40 / 40\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  #  Labels                                                 Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "  1  pore.back                                                               5\n",
      "  2  pore.domain_1                                                          25\n",
      "  3  pore.front                                                              5\n",
      "  4  pore.left                                                               5\n",
      "  5  pore.right                                                              5\n",
      "  6  pore.surface                                                           16\n",
      "  7  throat.domain_1                                                        40\n",
      "  8  throat.surface                                                         16\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f10eed9",
   "metadata": {},
   "source": [
    "Note the ``'throat.hydraulic_size_factors'`` array.  This is computed by a pore-scale model on the ``Demo`` network.  This computation is more rigorous in the following ways:\n",
    "\n",
    "1. The conductance of each half pore and that throat is considered.\n",
    "1. The throat length is computed carefully by accounting for the 'lens' between the intersection of the spherical pore bodies and the cylindrical throat.\n",
    "1. The net cross-sectional area of the pores are computed by integrating between the pore center and the pore-throat intersection point\n",
    "\n",
    "The conductance of each element in the conduit is returned as an *Nt-by-3* array, where columns 1 and 3 contain the hydraulic conductance of the half pore on either end of the throat, and the column 1 contains the throat conductance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "65e69a58",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[7.17224268e-17, 8.95163251e-19, 3.41740957e-17],\n",
       "       [3.41740957e-17, 9.39044534e-19, 8.24868939e-17],\n",
       "       [1.25531646e-16, 2.37360178e-18, 6.67669730e-17],\n",
       "       [6.67669730e-17, 2.40283561e-18, 1.29499002e-16],\n",
       "       [1.77366727e-16, 8.56740609e-18, 1.85484559e-16],\n",
       "       [2.00896764e-16, 1.22106603e-17, 3.02192896e-16],\n",
       "       [6.01031396e-16, 5.95606954e-17, 6.03870375e-16],\n",
       "       [6.09032656e-16, 6.33038224e-17, 6.44196320e-16],\n",
       "       [1.03665992e-16, 1.78551708e-18, 5.54807336e-17],\n",
       "       [5.54807336e-17, 1.79169742e-18, 1.04696234e-16],\n",
       "       [1.54582938e-16, 4.28175177e-18, 1.03819232e-16],\n",
       "       [1.03819232e-16, 3.85397995e-18, 1.10662774e-16],\n",
       "       [1.74675531e-16, 2.85535840e-18, 6.83639178e-17],\n",
       "       [6.83639178e-17, 2.51234835e-18, 1.35228170e-16],\n",
       "       [1.77395616e-16, 4.61695788e-18, 1.05679867e-16],\n",
       "       [1.05679867e-16, 4.28604303e-18, 1.46233002e-16],\n",
       "       [1.49269368e-16, 7.63835793e-18, 2.28705429e-16],\n",
       "       [1.05688678e-16, 1.35839458e-18, 4.33964312e-17],\n",
       "       [4.33964312e-17, 1.10265282e-18, 5.10502986e-17],\n",
       "       [6.76662146e-17, 2.20907752e-18, 9.62504448e-17],\n",
       "       [2.00160225e-16, 8.80767939e-18, 1.77366727e-16],\n",
       "       [3.41740957e-17, 8.74561874e-19, 6.63512743e-17],\n",
       "       [3.60328718e-16, 2.82338438e-17, 4.34693093e-16],\n",
       "       [6.67669730e-17, 2.58750204e-18, 1.52739862e-16],\n",
       "       [3.91422673e-16, 3.31532469e-17, 4.87753350e-16],\n",
       "       [1.77366727e-16, 9.03098345e-18, 2.13610710e-16],\n",
       "       [8.97176659e-17, 1.70496547e-18, 5.54807336e-17],\n",
       "       [3.89612755e-16, 2.19273313e-17, 3.02360428e-16],\n",
       "       [2.01412829e-16, 4.80596819e-18, 1.03819232e-16],\n",
       "       [2.38528795e-16, 6.37011185e-18, 1.23932048e-16],\n",
       "       [2.94518073e-16, 2.28948023e-17, 4.35296440e-16],\n",
       "       [5.54807336e-17, 1.54494686e-18, 5.98270528e-17],\n",
       "       [3.02360428e-16, 2.04032913e-17, 3.41440656e-16],\n",
       "       [1.03819232e-16, 3.79486014e-18, 1.04478685e-16],\n",
       "       [1.23932048e-16, 5.35260627e-18, 1.61465149e-16],\n",
       "       [2.85250597e-16, 8.61639662e-18, 1.49269368e-16],\n",
       "       [6.83639178e-17, 2.55471368e-18, 1.40619360e-16],\n",
       "       [1.01687811e-16, 1.33664490e-18, 4.33964312e-17],\n",
       "       [7.96485950e-17, 2.10315731e-18, 6.76662146e-17],\n",
       "       [2.00119660e-16, 8.65180085e-18, 1.74681839e-16]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pn['throat.hydraulic_size_factors']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29f6fccb",
   "metadata": {},
   "source": [
    "This data is called the ``size factor`` because it is purely the geometrical information required for the computation of the hydraulic conductance.  So the Hagan-Poisseiulle equation for each element is written as:\n",
    "\n",
    "$$ Q = \\frac{F_h}{\\mu} \\Delta P = g_h \\Delta P$$\n",
    "\n",
    "Note that both the $\\frac{\\pi R^4}{8}$ term and $L$ have been rolled into the $F_h$ value.\n",
    "\n",
    "The total conductance of the pore-throat-pore conduit can be found as the sum of three resistors in series.  Since we have conductance values, we add the inverses, and invert again.  The full expression for the hydraulic conductance between pores i and j, through throat k, is:\n",
    "\n",
    "\n",
    "$$ Q = \\bigg( \\frac{\\mu}{F_{h, i}} + \\frac{\\mu}{F_{h, k}} + \\frac{\\mu}{F_{h, j}} \\bigg) ^ {-1} \\Delta P $$\n",
    "\n",
    "This can be computed by hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0b26948e",
   "metadata": {},
   "outputs": [],
   "source": [
    "F_h = water['throat.hydraulic_size_factors']\n",
    "water['throat.hydraulic_conductance'] = (mu * (1/F_h).sum(axis=1))**(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "fc49009f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9.64890861e-16, 1.01200775e-15, 2.52019344e-15, 2.55102875e-15,\n",
       "       8.76379447e-15, 1.24146226e-14, 5.56745682e-14, 5.89529400e-14,\n",
       "       1.90491669e-15, 1.91150861e-15, 4.48459671e-15, 4.02523528e-15,\n",
       "       3.02123137e-15, 2.66531163e-15, 4.83218689e-15, 4.48520822e-15,\n",
       "       7.88493570e-15, 1.45652099e-15, 1.17908308e-15, 2.34297724e-15,\n",
       "       9.01642545e-15, 9.42596730e-16, 2.76479432e-14, 2.74408781e-15,\n",
       "       3.22015299e-14, 9.24896606e-15, 1.81841007e-15, 2.17482514e-14,\n",
       "       5.02794728e-15, 6.61517491e-15, 2.26770426e-14, 1.64158955e-15,\n",
       "       2.02647323e-14, 3.96006729e-15, 5.56764422e-15, 8.86707390e-15,\n",
       "       2.70972049e-15, 1.43348763e-15, 2.22665273e-15, 8.86413983e-15])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "water['throat.hydraulic_conductance']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3ef96480",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                       \r"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([153381.08855867, 160388.98970566, 173533.16165927, 166629.27675144,\n",
       "       153163.58781883])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sf = op.algorithms.StokesFlow(network=pn, phase=water)\n",
    "sf.set_value_BC(pores=pn.pores('left'), values=100_000)\n",
    "sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)\n",
    "soln = sf.run()\n",
    "sf['pore.pressure'][pn.pores('right')]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d59ec79",
   "metadata": {},
   "source": [
    "As can be seen the numbers are about the same as with the simple case, but should be somewhat more correct.  In fact, these above pressures are a bit higher, which is because the total conductance of the conduit is lower due to the inclusion of the pore body lengths into the total length, compared to above where only the throat length was included."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4535fccd",
   "metadata": {},
   "source": [
    "### Pore-Scale Model Method\n",
    "\n",
    "Instead of computing the hydraulic conductance manually as done above, there is a pore-scale model available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9891436d",
   "metadata": {},
   "outputs": [],
   "source": [
    "water.add_model(propname='throat.hydraulic_conductance',\n",
    "                model=op.models.physics.hydraulic_conductance.generic_hydraulic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fd34b943",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                       \r"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([153381.08855867, 160388.98970566, 173533.16165927, 166629.27675144,\n",
       "       153163.58781883])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sf = op.algorithms.StokesFlow(network=pn, phase=water)\n",
    "sf.set_value_BC(pores=pn.pores('left'), values=100_000)\n",
    "sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)\n",
    "soln = sf.run()\n",
    "sf['pore.pressure'][pn.pores('right')]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56ab903f",
   "metadata": {},
   "source": [
    "Which gives exactly the same result, without have to manually deal with the conductances-in-series calculation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
